MARKOV PROCESSES CONDITIONED ON THEIR LOCATION AT LARGE EXPONENTIAL TIMES.

Suppose that (Xt ) t ≥0 is a one-dimensional Brownian motion with negative drift -µ. It is possible to make sense of conditioning this process to be in the state 0 at an independent exponential random time and if we kill the conditioned process at the exponential time the resulting process is Markov. If we let the rate parameter of the random time go to 0, then the limit of the killed Markov process evolves like X conditioned to hit 0, after which time it behaves as X killed at the last time X visits 0. Equivalently, the limit process has the dynamics of the killed "bang-bang" Brownian motion that evolves like Brownian motion with positive drift +µ when it is negative, like Brownian motion with negative drift -µ when it is positive, and is killed according to the local time spent at 0. An extension of this result holds in great generality for a Borel right process conditioned to be in some state a at an exponential random time, at which time it is killed. Our proofs involve understanding the Campbell measures associated with local times, the use of excursion theory, and the development of a suitable analogue of the "bang-bang" construction for a general Markov process. As examples, we consider the special case when the transient Borel right process is a one-dimensional diffusion. Characterizing the limiting conditioned and killed process via its infinitesimal generator leads to an investigation of the h-transforms of transient one-dimensional diffusion processes that goes beyond what is known and is of independent interest.

RECOVERING A TREE FROM THE LENGTHS OF SUBTREES SPANNED BY A RANDOMLY CHOSEN SEQUENCE OF LEAVES.

Given an edge-weighted tree T with n leaves, sample the leaves uniformly at random without replacement and let Wk , 2 ≤ k ≤ n, be the length of the subtree spanned by the first k leaves. We consider the question, "Can T be identified (up to isomorphism) by the joint probability distribution of the random vector (W2, …, Wn )?" We show that if T is known a priori to belong to one of various families of edge-weighted trees, then the answer is, "Yes." These families include the edge-weighted trees with edge-weights in general position, the ultrametric edge-weighted trees, and certain families with equal weights on all edges such as (k + 1)-valent and rooted k-ary trees for k ≥ 2 and caterpillars.

DOOB-MARTIN COMPACTIFICATION OF A MARKOV CHAIN FOR GROWING RANDOM WORDS SEQUENTIALLY.

We consider a Markov chain that iteratively generates a sequence of random finite words in such a way that the nth word is uniformly distributed over the set of words of length 2n in which n letters are a and n letters are b: at each step an a and a b are shuffled in uniformly at random among the letters of the current word. We obtain a concrete characterization of the Doob-Martin boundary of this Markov chain and thereby delineate all the ways in which the Markov chain can be conditioned to behave at large times. Writing N(u) for the number of letters a (equivalently, b) in the finite word u, we show that a sequence (un ) n∈ℕ of finite words converges to a point in the boundary if, for an arbitrary word ν, there is convergence as n tends to infinity of the probability that the selection of N(ν) letters a and N(ν) letters b uniformly at random from un and maintaining their relative order results in ν. We exhibit a bijective correspondence between the points in the boundary and ergodic random total orders on the set {a1, b1, a2, b2, …} that have distributions which are separately invariant under finite permutations of the indices of the a's and those of the b's. We establish a further bijective correspondence between the set of such random total orders and the set of pairs (µ, ν) of diffuse probability measures on [0,1] such that ½(µ + ν) is Lebesgue measure: the restriction of the random total order to {a1, b1,…, an, bn } is obtained by taking X1,…, Xn (resp. Y1,… ,Yn ) i.i.d. with common distribution µ (resp. ν), letting (Z1,…, Z2n) be {X1, Y1,…, Xn , Yn } in increasing order, and declaring that the kth smallest element in the restricted total order is ai (resp. bj ) if Zk = Xi (resp. Zk = Yj ).

THE SEMIGROUP OF METRIC MEASURE SPACES AND ITS INFINITELY DIVISIBLE PROBABILITY MEASURES.

A metric measure space is a complete, separable metric space equipped with a probability measure that has full support. Two such spaces are equivalent if they are isometric as metric spaces via an isometry that maps the probability measure on the first space to the probability measure on the second. The resulting set of equivalence classes can be metrized with the Gromov-Prohorov metric of Greven, Pfaffelhuber and Winter. We consider the natural binary operation ⊞ on this space that takes two metric measure spaces and forms their Cartesian product equipped with the sum of the two metrics and the product of the two probability measures. We show that the metric measure spaces equipped with this operation form a cancellative, commutative, Polish semigroup with a translation invariant metric. There is an explicit family of continuous semicharacters that is extremely useful for, inter alia, establishing that there are no infinitely divisible elements and that each element has a unique factorization into prime elements. We investigate the interaction between the semigroup structure and the natural action of the positive real numbers on this space that arises from scaling the metric. For example, we show that for any given positive real numbers a, b, c the trivial space is the only space that satisfies a ⊞ b = c . We establish that there is no analogue of the law of large numbers: if X1, X2, … is an identically distributed independent sequence of random spaces, then no subsequence of [Formula: see text] converges in distribution unless each Xk is almost surely equal to the trivial space. We characterize the infinitely divisible probability measures and the Lévy processes on this semigroup, characterize the stable probability measures and establish a counterpart of the LePage representation for the latter class.

Bayesian Inference of Natural Selection from Allele Frequency Time Series.

The advent of accessible ancient DNA technology now allows the direct ascertainment of allele frequencies in ancestral populations, thereby enabling the use of allele frequency time series to detect and estimate natural selection. Such direct observations of allele frequency dynamics are expected to be more powerful than inferences made using patterns of linked neutral variation obtained from modern individuals. We developed a Bayesian method to make use of allele frequency time series data and infer the parameters of general diploid selection, along with allele age, in nonequilibrium populations. We introduce a novel path augmentation approach, in which we use Markov chain Monte Carlo to integrate over the space of allele frequency trajectories consistent with the observed data. Using simulations, we show that this approach has good power to estimate selection coefficients and allele age. Moreover, when applying our approach to data on horse coat color, we find that ignoring a relevant demographic history can significantly bias the results of inference. Our approach is made available in a C++ software package.

Protected polymorphisms and evolutionary stability of patch-selection strategies in stochastic environments.

We consider a population living in a patchy environment that varies stochastically in space and time. The population is composed of two morphs (that is, individuals of the same species with different genotypes). In terms of survival and reproductive success, the associated phenotypes differ only in their habitat selection strategies. We compute invasion rates corresponding to the rates at which the abundance of an initially rare morph increases in the presence of the other morph established at equilibrium. If both morphs have positive invasion rates when rare, then there is an equilibrium distribution such that the two morphs coexist; that is, there is a protected polymorphism for habitat selection. Alternatively, if one morph has a negative invasion rate when rare, then it is asymptotically displaced by the other morph under all initial conditions where both morphs are present. We refine the characterization of an evolutionary stable strategy for habitat selection from Schreiber (Am Nat 180:17-34, 2012) in a mathematically rigorous manner. We provide a necessary and sufficient condition for the existence of an ESS that uses all patches and determine when using a single patch is an ESS. We also provide an explicit formula for the ESS when there are two habitat types. We show that adding environmental stochasticity results in an ESS that, when compared to the ESS for the corresponding model without stochasticity, spends less time in patches with larger carrying capacities and possibly makes use of sink patches, thereby practicing a spatial form of bet hedging.

Evolutionary shaping of demographic schedules.

Evolutionary processes of natural selection may be expected to leave their mark on age patterns of survival and reproduction. Demographic theory includes three main strands--mutation accumulation, stochastic vitality, and optimal life histories. This paper reviews the three strands and, concentrating on mutation accumulation, extends a mathematical result with broad implications concerning the effect of interactions between small age-specific effects of deleterious mutant alleles. Empirical data from genomic sequencing along with prospects for combining strands of theory hold hope for future progress.

Analysis and rejection sampling of Wright-Fisher diffusion bridges.

We investigate the properties of a Wright-Fisher diffusion process starting at frequency x at time 0 and conditioned to be at frequency y at time T. Such a process is called a bridge. Bridges arise naturally in the analysis of selection acting on standing variation and in the inference of selection from allele frequency time series. We establish a number of results about the distribution of neutral Wright-Fisher bridges and develop a novel rejection-sampling scheme for bridges under selection that we use to study their behavior.

The age-specific force of natural selection and biodemographic walls of death.

W. D. Hamilton's celebrated formula for the age-specific force of natural selection furnishes predictions for senescent mortality due to mutation accumulation, at the price of reliance on a linear approximation. Applying to Hamilton's setting the full nonlinear demographic model for mutation accumulation recently developed by Evans, Steinsaltz, and Wachter, we find surprising differences. Nonlinear interactions cause the collapse of Hamilton-style predictions in the most commonly studied case, refine predictions in other cases, and allow walls of death at ages before the end of reproduction. Haldane's principle for genetic load has an exact but unfamiliar generalization.

Edge principal components and squash clustering: using the special structure of phylogenetic placement data for sample comparison.

Principal components analysis (PCA) and hierarchical clustering are two of the most heavily used techniques for analyzing the differences between nucleic acid sequence samples taken from a given environment. They have led to many insights regarding the structure of microbial communities. We have developed two new complementary methods that leverage how this microbial community data sits on a phylogenetic tree. Edge principal components analysis enables the detection of important differences between samples that contain closely related taxa. Each principal component axis is a collection of signed weights on the edges of the phylogenetic tree, and these weights are easily visualized by a suitable thickening and coloring of the edges. Squash clustering outputs a (rooted) clustering tree in which each internal node corresponds to an appropriate "average" of the original samples at the leaves below the node. Moreover, the length of an edge is a suitably defined distance between the averaged samples associated with the two incident nodes, rather than the less interpretable average of distances produced by UPGMA, the most widely used hierarchical clustering method in this context. We present these methods and illustrate their use with data from the human microbiome.

Stochastic population growth in spatially heterogeneous environments.

Classical ecological theory predicts that environmental stochasticity increases extinction risk by reducing the average per-capita growth rate of populations. For sedentary populations in a spatially homogeneous yet temporally variable environment, a simple model of population growth is a stochastic differential equation dZ(t) = µZ(t)dt + σZ(t)dW(t), t ≥ 0, where the conditional law of Z(t+Δt)-Z(t) given Z(t) = z has mean and variance approximately z µΔt and z²σ²Δt when the time increment Δt is small. The long-term stochastic growth rate lim(t→∞) t⁻¹ log Z(t) for such a population equals µ − σ²/2 . Most populations, however, experience spatial as well as temporal variability. To understand the interactive effects of environmental stochasticity, spatial heterogeneity, and dispersal on population growth, we study an analogous model X(t) = (X¹(t) , . . . , X(n)(t)), t ≥ 0, for the population abundances in n patches: the conditional law of X(t+Δt) given X(t) = x is such that the conditional mean of X(i)(t+Δt) − X(i)(t) is approximately [x(i)µ(i) + Σ(j) (x(j) D(ji) − x(i) D(i j) )]Δt where µ(i) is the per capita growth rate in the ith patch and D(ij) is the dispersal rate from the ith patch to the jth patch, and the conditional covariance of X(i)(t+Δt)− X(i)(t) and X(j)(t+Δt) − X(j)(t) is approximately x(i)x(j)σ(ij)Δt for some covariance matrix Σ = (σ(ij)). We show for such a spatially extended population that if S(t) = X¹(t)+· · ·+ X(n)(t) denotes the total population abundance, then Y(t) = X(t)/S(t), the vector of patch proportions, converges in law to a random vector Y(∞) as t → ∞, and the stochastic growth rate lim(t→∞) t⁻¹ log S(t) equals the space-time average per-capita growth rate Σ(i)µ(i)E[Y(i)(∞)] experienced by the population minus half of the space-time average temporal variation E[Σ(i,j) σ(i j)Y(i)(∞) Y(j)(∞)] experienced by the population. Using this characterization of the stochastic growth rate, we derive an explicit expression for the stochastic growth rate for populations living in two patches, determine which choices of the dispersal matrix D produce the maximal stochastic growth rate for a freely dispersing population, derive an analytic approximation of the stochastic growth rate for dispersal limited populations, and use group theoretic techniques to approximate the stochastic growth rate for populations living in multi-scale landscapes (e.g. insects on plants in meadows on islands). Our results provide fundamental insights into "ideal free" movement in the face of uncertainty, the persistence of coupled sink populations, the evolution of dispersal rates, and the single large or several small (SLOSS) debate in conservation biology. For example, our analysis implies that even in the absence of density-dependent feedbacks, ideal-free dispersers occupy multiple patches in spatially heterogeneous environments provided environmental fluctuations are sufficiently strong and sufficiently weakly correlated across space. In contrast, for diffusively dispersing populations living in similar environments, intermediate dispersal rates maximize their stochastic growth rate.

Estimating allele age and selection coefficient from time-serial data.

Recent advances in sequencing technologies have made available an ever-increasing amount of ancient genomic data. In particular, it is now possible to target specific single nucleotide polymorphisms in several samples at different time points. Such time-series data are also available in the context of experimental or viral evolution. Time-series data should allow for a more precise inference of population genetic parameters and to test hypotheses about the recent action of natural selection. In this manuscript, we develop a likelihood method to jointly estimate the selection coefficient and the age of an allele from time-serial data. Our method can be used for allele frequencies sampled from a single diallelic locus. The transition probabilities are calculated by approximating the standard diffusion equation of the Wright-Fisher model with a one-step process. We show that our method produces unbiased estimates. The accuracy of the method is tested via simulations. Finally, the utility of the method is illustrated with an application to several loci encoding coat color in horses, a pattern that has previously been linked with domestication. Importantly, given our ability to estimate the age of the allele, it is possible to gain traction on the important problem of distinguishing selection on new mutations from selection on standing variation. In this coat color example for instance, we estimate the age of this allele, which is found to predate domestication.

The phylogenetic Kantorovich-Rubinstein metric for environmental sequence samples.

It is now common to survey microbial communities by sequencing nucleic acid material extracted in bulk from a given environment. Comparative methods are needed that indicate the extent to which two communities differ given data sets of this type. UniFrac, which gives a somewhat ad hoc phylogenetics-based distance between two communities, is one of the most commonly used tools for these analyses. We provide a foundation for such methods by establishing that, if we equate a metagenomic sample with its empirical distribution on a reference phylogenetic tree, then the weighted UniFrac distance between two samples is just the classical Kantorovich-Rubinstein, or earth mover's, distance between the corresponding empirical distributions. We demonstrate that this Kantorovich-Rubinstein distance and extensions incorporating uncertainty in the sample locations can be written as a readily computable integral over the tree, we develop L(p) Zolotarev-type generalizations of the metric, and we show how the p-value of the resulting natural permutation test of the null hypothesis 'no difference between two communities' can be approximated by using a Gaussian process functional. We relate the L(2)-case to an analysis-of-variance type of decomposition, finding that the distribution of its associated Gaussian functional is that of a computable linear combination of independent [Formula: see text] random variables.

Ubiquity of synonymity: almost all large binary trees are not uniquely identified by their spectra or their immanantal polynomials.

BACKGROUND: There are several common ways to encode a tree as a matrix, such as the adjacency matrix, the Laplacian matrix (that is, the infinitesimal generator of the natural random walk), and the matrix of pairwise distances between leaves. Such representations involve a specific labeling of the vertices or at least the leaves, and so it is natural to attempt to identify trees by some feature of the associated matrices that is invariant under relabeling. An obvious candidate is the spectrum of eigenvalues (or, equivalently, the characteristic polynomial). RESULTS: We show for any of these choices of matrix that the fraction of binary trees with a unique spectrum goes to zero as the number of leaves goes to infinity. We investigate the rate of convergence of the above fraction to zero using numerical methods. For the adjacency and Laplacian matrices, we show that the a priori more informative immanantal polynomials have no greater power to distinguish between trees. CONCLUSION: Our results show that a generic large binary tree is highly unlikely to be identified uniquely by common spectral invariants.

Transcriptional regulation: effects of promoter proximal pausing on speed, synchrony and reliability.

Recent whole genome polymerase binding assays in the Drosophila embryo have shown that a substantial proportion of uninduced genes have pre-assembled RNA polymerase-II transcription initiation complex (PIC) bound to their promoters. These constitute a subset of promoter proximally paused genes for which mRNA elongation instead of promoter access is regulated. This difference can be described as a rearrangement of the regulatory topology to control the downstream transcriptional process of elongation rather than the upstream transcriptional initiation event. It has been shown experimentally that genes with the former mode of regulation tend to induce faster and more synchronously, and that promoter-proximal pausing is observed mainly in metazoans, in accord with a posited impact on synchrony. However, it has not been shown whether or not it is the change in the regulated step per se that is causal. We investigate this question by proposing and analyzing a continuous-time Markov chain model of PIC assembly regulated at one of two steps: initial polymerase association with DNA, or release from a paused, transcribing state. Our analysis demonstrates that, over a wide range of physical parameters, increased speed and synchrony are functional consequences of elongation control. Further, we make new predictions about the effect of elongation regulation on the consistent control of total transcript number between cells. We also identify which elements in the transcription induction pathway are most sensitive to molecular noise and thus possibly the most evolutionarily constrained. Our methods produce symbolic expressions for quantities of interest with reasonable computational effort and they can be used to explore the interplay between interaction topology and molecular noise in a broader class of biochemical networks. We provide general-purpose code implementing these methods.

Shape-based peak identification for ChIP-Seq.

BACKGROUND: The identification of binding targets for proteins using ChIP-Seq has gained popularity as an alternative to ChIP-chip. Sequencing can, in principle, eliminate artifacts associated with microarrays, and cheap sequencing offers the ability to sequence deeply and obtain a comprehensive survey of binding. A number of algorithms have been developed to call "peaks" representing bound regions from mapped reads. Most current algorithms incorporate multiple heuristics, and despite much work it remains difficult to accurately determine individual peaks corresponding to distinct binding events. RESULTS: Our method for identifying statistically significant peaks from read coverage is inspired by the notion of persistence in topological data analysis and provides a non-parametric approach that is statistically sound and robust to noise in experiments. Specifically, our method reduces the peak calling problem to the study of tree-based statistics derived from the data. We validate our approach using previously published data and show that it can discover previously missed regions. CONCLUSIONS: The difficulty in accurately calling peaks for ChIP-Seq data is partly due to the difficulty in defining peaks, and we demonstrate a novel method that improves on the accuracy of previous methods in resolving peaks. Our introduction of a robust statistical test based on ideas from topological data analysis is also novel. Our methods are implemented in a program called T-PIC (Tree shape Peak Identification for ChIP-Seq) is available at http://bio.math.berkeley.edu/tpic/.

Coverage statistics for sequence census methods.

BACKGROUND: We study the statistical properties of fragment coverage in genome sequencing experiments. In an extension of the classic Lander-Waterman model, we consider the effect of the length distribution of fragments. We also introduce a coding of the shape of the coverage depth function as a tree and explain how this can be used to detect regions with anomalous coverage. This modeling perspective is especially germane to current high-throughput sequencing experiments, where both sample preparation protocols and sequencing technology particulars can affect fragment length distributions. RESULTS: Under the mild assumptions that fragment start sites are Poisson distributed and successive fragment lengths are independent and identically distributed, we observe that, regardless of fragment length distribution, the fragments produced in a sequencing experiment can be viewed as resulting from a two-dimensional spatial Poisson process. We then study the successive jumps of the coverage function, and show that they can be encoded as a random tree that is approximately a Galton-Watson tree with generation-dependent geometric offspring distributions whose parameters can be computed. CONCLUSIONS: We extend standard analyses of shotgun sequencing that focus on coverage statistics at individual sites, and provide a null model for detecting deviations from random coverage in high-throughput sequence census based experiments. Our approach leads to explicit determinations of the null distributions of certain test statistics, while for others it greatly simplifies the approximation of their null distributions by simulation. Our focus on fragments also leads to a new approach to visualizing sequencing data that is of independent interest.

Confidence intervals for negative binomial random variables of high dispersion.

We consider the problem of constructing confidence intervals for the mean of a Negative Binomial random variable based upon sampled data. When the sample size is large, it is a common practice to rely upon a Normal distribution approximation to construct these intervals. However, we demonstrate that the sample mean of highly dispersed Negative Binomials exhibits a slow convergence in distribution to the Normal as a function of the sample size. As a result, standard techniques (such as the Normal approximation and bootstrap) will construct confidence intervals for the mean that are typically too narrow and significantly undercover at small sample sizes or high dispersions. To address this problem, we propose techniques based upon Bernstein's inequality or the Gamma and Chi Square distributions as alternatives to the standard methods. We investigate the impact of imposing a heuristic assumption of boundedness on the data as a means of improving the Bernstein method. Furthermore, we propose a ratio statistic relating the Negative Binomial's parameters that can be used to ascertain the applicability of the Chi Square method and to provide guidelines on evaluating the length of all proposed methods. We compare the proposed methods to the standard techniques in a variety of simulation experiments and consider data arising in the serial analysis of gene expression and traffic flow in a communications network.

Vital Rates from the Action of Mutation Accumulation.

New models for evolutionary processes of mutation accumulation allow hypotheses about the age-specificity of mutational effects to be translated into predictions of heterogeneous population hazard functions. We apply these models to questions in the biodemography of longevity, including proposed explanations of Gompertz hazards and mortality plateaus.

To what extent does genealogical ancestry imply genetic ancestry?

Recent statistical and computational analyses have shown that a genealogical most recent common ancestor (MRCA) may have lived in the recent past [Chang, J.T., 1999. Recent common ancestors of all present-day individuals. Adv. Appl. Probab. 31, 1002-1026. 1027-1038; Rohde, D.L.T., Olson, S., Chang, J.T., 2004. Modelling the recent common ancestry of all living humans. Nature 431, 562-566]. However, coalescent-based approaches show that genetic most recent common ancestors for a given non-recombining locus are typically much more ancient [Kingman, J.F.C., 1982a. The coalescent. Stochastic Process Appl. 13, 235-248; Kingman, J.F.C., 1982b. On the genealogy of large populations. J. Appl. Probab. 19A, 27-43]. It is not immediately clear how these two perspectives interact. This paper investigates relationships between the number of descendant alleles of an ancestor allele and the number of genealogical descendants of the individual who possessed that allele for a simple diploid genetic model extending the genealogical model of [Chang, J.T., 1999. Recent common ancestors of all present-day individuals. Adv. Appl. Probab. 31, 1002-1026. 1027-1038].